{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "78da6786",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dbfbf93a",
   "metadata": {},
   "outputs": [],
   "source": [
    "def refinement_point(newKnots,weightedPts,knotVec,p,oldKnotVec):\n",
    "    for knot in range(len(newKnots)):  # loop over all the new knots\n",
    "        newControlPts=np.zeros(len(weightedPts)+1,3)\n",
    "        knot_bar     =newKnots[knot]\n",
    "        kval         =np.where(knotVec>knot_bar)[0]-1\n",
    "        for i in range(len(newControlPts)):\n",
    "            if i<= (kval - p):\n",
    "                alpha=1\n",
    "            elif i>=(kval-p+1) and i<= kval:\n",
    "                alpha= (newKnots[knot]-oldKnotVec[i])/(oldKnotVec[i+p]-oldKnotVec[i])\n",
    "            else:\n",
    "                alpha=0\n",
    "            newPoints=np.zeros(1,3)\n",
    "            if i != 1:\n",
    "                newPoints=(1-alpha)*weightedPts[i-1,:]\n",
    "            if i != len(newControlPts):\n",
    "                newPoints=newPoints + alpha * weightedPts[i,:]\n",
    "            newControlPts[i,:] = newPoints\n",
    "        weightedPts = newControlPts\n",
    "        knotVec     = np.sort([knotVec, knot_bar])\n",
    "        oldKnotVec  = knotVec\n",
    "    controlPts =np.array([weightedPts[:,0]/weightedPts[:,2], weightedPts[:, 1]/weightedPts[:, 2], weightedPts[:, 2]])\n",
    "    return controlPts\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ed5aa5c2",
   "metadata": {},
   "outputs": [],
   "source": [
    "def h_refinement( p,q,refinementu,refinementv,knotVec_u,knotVec_v,controlPts):\n",
    "    uniqueKnots_u = np.unique(knotVec_u)\n",
    "    uniqueKnots_v = np.unique(knotVec_v)\n",
    "    weightedPts = np.array([controlPts[:,0]*controlPts[:,2], controlPts[:,2]*controlPts[:,2], controlPts[:,2]])\n",
    "    \n",
    "    oldKnotVec_u=knotVec_u\n",
    "    oldKnotVec_v=knotVec_v\n",
    "    newKnots_u  = np.zeros((len(uniqueKnots_u)-1),refinementu-1) \n",
    "    newKnots_v  = np.zeros((len(uniqueKnots_v)-1),refinementv-1)\n",
    "    for i in range(len(newKnots_u)):\n",
    "        distributed_u  = np.linspace(uniqueKnots_u[i],uniqueKnots_u[i+1], refinementu +1)\n",
    "        newKnots_u[i,:]=distributed_u[2:end-1]\n",
    "    for i in range(len(newKnots_v)):\n",
    "        distributed_v  =np.linspace(uniqueKnots_v[i],uniqueKnots_v[i+1], refinementv+1)\n",
    "        newKnots_v[i,:]=distributed_v[1:end-2]\n",
    "    mu,nu=newKnots_u.shape\n",
    "    newKnots_u = np.sort(newKnots_u.reshape(1, mu*nu))\n",
    "    newKnotsu= np.sort([newKnots_u, knotVec_u])\n",
    "    mv,nv=newKnots_v.shape\n",
    "    newKnots_v=np.sort(newKnots_v.reshape(1,mv*nv))\n",
    "    newKnotsv=np.sort([newKnots_v, knotVec_v])\n",
    "    \n",
    "    # and now redefine the control points\n",
    "    numBasisFns_u=len(knotVec_u)-1-p\n",
    "    numBasisFns_v=len(knotVec_v)-1-q\n",
    "    \n",
    "    tcontrolPts  = refinement_point(newKnots_u, weightedPts[:numBasisFns_u,:],knotVec_u,p,oldKnotVec_u)\n",
    "    ttcontrolPts = refinement_point(newKnots_v,\n",
    "            weightedPts[:numBasisFns_u:(numBasisFns_v-1)*numBasisFns_u+1,:],\n",
    "            knotVec_v,q,oldKnotVec_v)\n",
    "    trancontrolPts=np.zeros(tcontrolPts.shape[0], ttcontrolPts.shape[0], 3)\n",
    "    for i in range(numBasisFns_u):\n",
    "        ttcontrolPts = refinement_point(newKnots_v, \n",
    "            weightedPts[i:numBasisFns_u:(numBasisFns_v-1)*numBasisFns_u+i,:], knotVec_v, q, oldKnotVec_v)\n",
    "        trancontrolPts[i,:,:]=ttcontrolPts\n",
    "    \n",
    "    #对v方向进行细化\n",
    "    mm=np.zeros(numBasisFns_u,3)\n",
    "    for j in range(ttcontrolPts.shape[0]):\n",
    "        for jj in range(numBasisFns_u):\n",
    "            mm[jj,:]=trancontrolPts[jj,j,:]\n",
    "        mm =np.array([mm[:,1]*mm[:,3], mm[:,2]*mm[:,3], mm[:,3]])\n",
    "        tcontrolPts = refinement_point(newKnots_u,\n",
    "            mm,knotVec_u,p,oldKnotVec_u)\n",
    "        trancontrolPts[:,j,:] = tcontrolPts\n",
    "    controlPts=np.zeros(trancontrolPts.shape[0]*trancontrolPts.shape[1],3)\n",
    "    for i in range(trancontrolPts.shape[0]):\n",
    "        for j in range(trancontrolPts.shape[1]):\n",
    "            controlPts[i+(j-1)*(trancontrolPts.shape[0]),:]=trancontrolPts[i,j,:]\n",
    "    return controlPts,newKnotsu, newKnotsv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "81773701",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# 初始模型参数 \n",
    "\n",
    "# l形梁\n",
    "knotVecu = [0,0, 0, 0.5, 0.5, 1, 1, 1]\n",
    "p=2\n",
    "knotVecv = [0, 0, 0, 1, 1, 1]\n",
    "q=2\n",
    "# knotVec_u = [0 0 0 0.5 0.5 1 1 1]\n",
    "# p=2\n",
    "# knotVec_v = [0 0 0 1 1 1]\n",
    "# q=2\n",
    "\n",
    "controlPts = np.array([\n",
    "    [-0.12,  0.18, 1],\n",
    "    [-0.12,  0.00,  1],\n",
    "    [-0.12, -0.12, 1],\n",
    "    [0,-0.12, 1],\n",
    "    [0.18, -0.12, 1],\n",
    "    [-0.06,0.18,1],\n",
    "    [-0.06,0,1],\n",
    "    [-0.06,-0.06,1],\n",
    "    [0,-0.06,1],\n",
    "    [0.18,-0.06,1],\n",
    "    [0,0.18,1],\n",
    "    [0,0.09,1],\n",
    "    [0,0,1],\n",
    "    [0.09,0,1],\n",
    "    [0.18,0,1],\n",
    "])\n",
    "\n",
    "# 细化次数 \n",
    "refinementu=5\n",
    "refinementv=5\n",
    "controlPts,knotVecu,knotVecv = h_refinement(p,q,refinementu,refinementv,knotVecu, knotVecv,controlPts ) # 单元细化"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8d5bfbdd",
   "metadata": {},
   "outputs": [],
   "source": [
    "def iso_numelement(knotVec,p):\n",
    "    uniqueKnots = np.unique(knotVec)\n",
    "    ne = len(uniqueKnots)-1   # number of elements\n",
    "    elRange = np.zeros(ne,2)        # initialise matrices\n",
    "    elConn  = np.zeros(ne,p+1) \n",
    "    elKnotIndices = np.zeros(ne,2)\n",
    "    tracConn      = np.zeros(ne,p+1)\n",
    "\n",
    "    # determine our element ranges and the corresponding knot indices\n",
    "    element =1\n",
    "    previousKnotVal=0\n",
    "    for i in range(len(knotVec)):\n",
    "        currentKnotVal=knotVec[i]\n",
    "        if knotVec[i] != previousKnotVal:\n",
    "            elRange[element,:]= np.array([previousKnotVal, currentKnotVal])\n",
    "            elKnotIndices[element,:]=np.array([i-1, i])\n",
    "            element=element+1\n",
    "        \n",
    "        previousKnotVal=currentKnotVal\n",
    "\n",
    "    numRepeatedKnots=0\n",
    "    for e in range(ne):\n",
    "        indices=np.arange(int(elKnotIndices[e,1]-p+1), int(elKnotIndices[e,1]))\n",
    "        previousKnotVals=knotVec[indices]\n",
    "        currentKnotVals=np.ones(1,p)*knotVec(elKnotIndices[e,1])\n",
    "        if previousKnotVals == currentKnotVals and np.sum(previousKnotVals != 0)>1:\n",
    "            numRepeatedKnots=numRepeatedKnots+1\n",
    "        \n",
    "        elConn[e,:]=np.arange((elKnotIndices(e,1)-p), elKnotIndices(e,1))\n",
    "        tracConn[e,:]=elConn[e,:]+numRepeatedKnots\n",
    "    \n",
    "    bsFnConn=elConn;       # 组成每个节点格子的两个节点索引\n",
    "    return bsFnConn,elRange\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dd5603ab",
   "metadata": {},
   "outputs": [],
   "source": [
    "def ISOnumtotal(knotVec_u,p,knotVec_v,q,controlPts):\n",
    "\n",
    "    # knotVec_u   =[0 0 0 0.5  1 1 1]\n",
    "    # knotVec_v   =[0 0 0  0.5 1 1 1]\n",
    "    # q=2\n",
    "    # p=2\n",
    "    # controlPts=[1 0 1;2.5 0 1;4 0 1;1 0.4142136 0.8536;2.5 0.75 1;4 4 1;\n",
    "    #     0.4142136 1 0.8536;0.75 2.5 1;4 4 1;0 1 1;0,2.5 1;0 4 1]\n",
    "\n",
    "    # refinementu=5\n",
    "    # refinementv=5\n",
    "    #  [ controlPts,knotVec_u,knotVec_v] = h_refinement( p,q,refinementu,refinementv,knotVec_u,...\n",
    "    #   knotVec_v,controlPts )\n",
    "\n",
    "    bsFnConn_u, _ =iso_numelement(knotVec_u,p)\n",
    "    bsFnConn_v, _ =iso_numelement(knotVec_v,q)\n",
    "    numBasisFns_u=len(knotVec_u)-1-p\n",
    "    numBasisFns_v=len(knotVec_v)-1-q\n",
    "    num=np.zeros(numBasisFns_v,numBasisFns_u)\n",
    "    for i in range(numBasisFns_v):\n",
    "        for j in range(numBasisFns_u):\n",
    "            num[i, j] = numBasisFns_u*(i-1)+j\n",
    "    numm=num\n",
    "    element=np.zeros(bsFnConn_u.shape[0]*bsFnConn_v.shape[0], bsFnConn_u.shape[1]*bsFnConn_v.shape[1])\n",
    "\n",
    "    for i in range(bsFnConn_u.shape[0]):\n",
    "        for j in range(bsFnConn_v.shape[0]):\n",
    "            element[(j-1)*bsFnConn_u.shape[0]+i, :]=np.sort(num[bsFnConn_v[j,:],bsFnConn_u[i,:]].reshape(1, bsFnConn_u.shape[1]*bsFnConn_v.shape[1]))\n",
    "\n",
    "\n",
    "    elenode=np.array([np.arange(numBasisFns_v*numBasisFns_u), controlPts])\n",
    "\n",
    "    return element,elenode,numm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ed2625b9",
   "metadata": {},
   "outputs": [],
   "source": [
    "def ISOnum(knotVecu,knotVecv,iel):\n",
    "    uniqueKnots_u = np.unique(knotVecu)\n",
    "    uniqueKnots_v = np.unique(knotVecv)\n",
    "\n",
    "\n",
    "    isonum=np.zeros(4,2)\n",
    "\n",
    "    numBasisFns_u=len(uniqueKnots_u)\n",
    "    numBasisFns_v=len(uniqueKnots_v)\n",
    "\n",
    "    eu=np.zeros(numBasisFns_u-1,2)\n",
    "    eu[:,0]=np.arange(1,numBasisFns_u-1).T\n",
    "    eu[:,1]=np.arange(2, numBasisFns_u).T\n",
    "    ev=np.zeros(numBasisFns_v-1,2)\n",
    "    ev[:,0] = np.arange(1, numBasisFns_v-1)\n",
    "    ev[:,1] = np.arange(2, numBasisFns_v)\n",
    "\n",
    "    # iu = mod(iel,numBasisFns_u-1)\n",
    "    iu = iel % numBasisFns_u-1\n",
    "    if iu==0:\n",
    "        iu=numBasisFns_u-1\n",
    "\n",
    "    iv=(iel-iu)/(numBasisFns_u-1)+1\n",
    "    isonum[0,0]=uniqueKnots_u(eu(iu,1))\n",
    "    isonum[0,1]=uniqueKnots_v(ev(iv,1))\n",
    "    isonum[1,0]=uniqueKnots_u(eu(iu,2))\n",
    "    isonum[1,1]=uniqueKnots_v(ev(iv,1))\n",
    "    isonum[2,0]=uniqueKnots_u(eu(iu,1))\n",
    "    isonum[2,1]=uniqueKnots_v(ev(iv,2))\n",
    "    isonum[3,0]=uniqueKnots_u(eu(iu,2))\n",
    "    isonum[3,1]=uniqueKnots_v(ev(iv,2))\n",
    "    return isonum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "814ae665",
   "metadata": {},
   "outputs": [],
   "source": [
    "def GaussPoint2(No_point_x, No_point_y):\n",
    "    # %--------------------------------------------------------------------------\n",
    "    # %  Purpose:\n",
    "    # %     determine the integration points and weighting coefficients\n",
    "    # %     of Gauss quadrature for two-dimensional integration\n",
    "    # %--------------------------------------------------------------------------\n",
    "    # %----------------------------------------------------------\n",
    "    # %  (1) determine the largest one between No_point_x and No_point_y\n",
    "    # %----------------------------------------------------------\n",
    "    if No_point_x > No_point_y:\n",
    "        ng = No_point_x\n",
    "    else:\n",
    "        ng = No_point_y\n",
    "    # %----------------------------------\n",
    "    # %  (2) initialization of coefficient vector\n",
    "    # %----------------------------------\n",
    "    point2 = np.zeros(ng, 2)\n",
    "    weight2 = np.zeros(ng, 2)\n",
    "    # %-----------------------------------------------\n",
    "    # %  (3) find corresponding integration points and weights\n",
    "    # %-----------------------------------------------\n",
    "    pointx, weightx = GaussPoint1(No_point_x)  # quadrature rule for x-axis\n",
    "    pointy, weighty = GaussPoint1(No_point_y)  # quadrature rule for y-axis\n",
    "    # %-------------------------------\n",
    "    # %  ��4) quadrature for two-dimension\n",
    "    # %-------------------------------\n",
    "    for ii in range(No_point_x):  # quadrature in x-axis\n",
    "        point2[ii, 0] = pointx(ii)\n",
    "        weight2[ii, 0] = weightx(ii)\n",
    "        for ij in range(No_point_y):  # quadrature in y-axis\n",
    "            point2[ij, 1] = pointy(ij)\n",
    "            weight2[ij, 1] = weighty(ij)\n",
    "    # %--------------------------------------------------------------------------\n",
    "    # %    The end\n",
    "    # %--------------------------------------------------------------------------\n",
    "\n",
    "    return point2, weight2\n",
    "\n",
    "\n",
    "def GaussPoint1(No_points):\n",
    "    # %--------------------------------------------------------------------------\n",
    "    # %  Purpose:\n",
    "    # %     determine the integration points and weighting coefficients of Gauss quadrature\n",
    "    # %  Variable Description:\n",
    "    # %     No_points - number of integration points\n",
    "    # %     point1  - vector containing integration points\n",
    "    # %     weight1 - vector containing weighting coefficients\n",
    "    # %-------------------------------------------------------------------\n",
    "    # %------------------------------------------\n",
    "    # %  (1) initialization of the vectors\n",
    "    # %------------------------------------------\n",
    "    point1 = np.zeros(No_points, 1)\n",
    "    weight1 = np.zeros(No_points, 1)\n",
    "\n",
    "    # %----------------------------------------------------------\n",
    "    # %  (2) find corresponding integration points and weights\n",
    "    # %----------------------------------------------------------\n",
    "    if No_points == 1:  # 1-point quadrature rule\n",
    "        point1[0] = 0.0\n",
    "        weight1[0] = 2.0\n",
    "    elif No_points == 2:  # 2-point quadrature rule\n",
    "        point1[0] = -0.577350269189626\n",
    "        point1[1] = -point1[0]\n",
    "\n",
    "        weight1[0] = 1.0\n",
    "        weight1[1] = weight1[0]\n",
    "    elif No_points == 3:  # 3-point quadrature rule\n",
    "        point1[0] = -0.774596669241483\n",
    "        point1[1] = 0.0\n",
    "        point1[2] = -point1[0]\n",
    "\n",
    "        weight1[0] = 0.555555555555556\n",
    "        weight1[1] = 0.888888888888889\n",
    "        weight1[2] = weight1[0]\n",
    "    elif No_points == 4:  # 4-point quadrature rule\n",
    "        point1[0] = -0.861136311594053\n",
    "        point1[1] = -0.339981043584856\n",
    "        point1[2] = -point1[1]\n",
    "        point1[3] = -point1[0]\n",
    "\n",
    "        weight1[0] = 0.347854845137454\n",
    "        weight1[1] = 0.652145154862546\n",
    "        weight1[2] = weight1[1]\n",
    "        weight1[3] = weight1[0]\n",
    "    elif No_points == 5:  # 5-point quadrature rule\n",
    "        point1[0] = -0.906179845938664\n",
    "        point1[1] = -0.538469310105683\n",
    "        point1[2] = 0.0\n",
    "        point1[3] = -point1[1]\n",
    "        point1[4] = -point1[0]\n",
    "\n",
    "        weight1[0] = 0.236926885056189\n",
    "        weight1[1] = 0.478628670499366\n",
    "        weight1[2] = 0.568888888888889\n",
    "        weight1[3] = weight1(2)\n",
    "        weight1[4] = weight1(1)\n",
    "    elif No_points == 6:  # 6-point quadrature rule\n",
    "        point1[0] = -0.932469514203152\n",
    "        point1[1] = -0.661209386466265\n",
    "        point1[2] = -0.238619186083197\n",
    "        point1[3] = -point1[2]\n",
    "        point1[4] = -point1[1]\n",
    "        point1[5] = -point1[0]\n",
    "\n",
    "        weight1[0] = 0.171324492379170\n",
    "        weight1[1] = 0.360761573048139\n",
    "        weight1[2] = 0.467913934572691\n",
    "        weight1[3] = weight1[2]\n",
    "        weight1[4] = weight1[1]\n",
    "        weight1[5] = weight1[0]\n",
    "    else:\n",
    "        print('number of integration points should be')\n",
    "        print('> 0 or < 7')\n",
    "    return point1, weight1\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8ee0d537",
   "metadata": {},
   "outputs": [],
   "source": [
    "def BsplineBasisHighOrder( lowerBasis, xi, knot, i, p ):\n",
    "    # the recursive function that gives us the basis functions for any p>0\n",
    "    eps = 2.2204e-16\n",
    "    if(lowerBasis[0] < eps):\n",
    "         basis=0\n",
    "    else:\n",
    "         basis=(xi-knot(i))/(knot(i+p)-knot(i)) * lowerBasis(1)\n",
    "    if(lowerBasis[1] > eps):\n",
    "        basis=basis + (knot(i+p+1) - xi)/(knot(i+1+p)-knot(i+1)) * lowerBasis(2)\n",
    "    return basis\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "99986494",
   "metadata": {},
   "outputs": [],
   "source": [
    "def BsplineBasis(i, p, xi, knot):\n",
    "    # % calculate the non-rational Bspline basis function\n",
    "\n",
    "    # % i=basis fn we want\n",
    "    # % p=order (0=const, 1=linear, 2=quadratic)\n",
    "    # % xi=parametric coordinate along b-spline\n",
    "    # % knot= the knot vector (eg. [0 0 0 1 2 3 4 4 5 5 5]) in parameter space\n",
    "    basis = np.zeros(p + 1, p + 1)  # this will contain the basis functions of everything of lower order\n",
    "    # % [ N_1,0 N_2,0 N_3,0 ..\n",
    "    # %   N_2,0 N_2,1 ...\n",
    "    # %   ETC.\n",
    "\n",
    "    # % first calculate the basis fns for order 0\n",
    "\n",
    "    # % The first if statement is the conventional interval for the zeroth order\n",
    "    # % B-spline functions. But if we want the curve to be interpolatory at the\n",
    "    # % end of the parameter space, we need to put another check in. Basically,\n",
    "    # % we want the zeroth order basis function for the max non-zero interval to\n",
    "    # % be equal to one at the end of the interval.\n",
    "    for c in range(p):\n",
    "        if (xi >= knot(i + c) and xi < knot(i + c + 1)):  # do we lie within the interval xi_i<=xi<xi_(i+!)?\n",
    "            basis[0, c] = 1\n",
    "        elif xi == max(knot) and knot(i + c + 1) == max(knot) and knot(i + c) != knot(i + c + 1):\n",
    "            basis[0, c] = 1\n",
    "\n",
    "        if (p == 0):  # if the basis order is zero we have the fn directly\n",
    "            N = basis[0, c]\n",
    "            dN = 0\n",
    "            return N, dN\n",
    "    # and now for the higher order fns\n",
    "    for b in range(1, p + 1):\n",
    "        start = p + 1 - b + 1\n",
    "        finish = 1\n",
    "        for c in range(finish, start):\n",
    "            basis[b, c] = BsplineBasisHighOrder(basis[b - 1, c:(c + 1)], xi, knot, i + c - 1, b - 1)\n",
    "\n",
    "    N = basis[p, 0]\n",
    "\n",
    "    temp = p / (knot[i + p] - knot[i]) * basis[p, 0]\n",
    "    if np.isnan(temp):\n",
    "        temp = 0\n",
    "    temp2 = p / (knot[i + p] - knot[i]) * basis[p, 1]\n",
    "    if np.isnan(temp2):\n",
    "        temp2 = 0\n",
    "    dN = temp - temp2\n",
    "\n",
    "    return N, dN\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "37867434",
   "metadata": {},
   "outputs": [],
   "source": [
    "def NURBSbasis_matlab(i, j, p, q, xi, yi, knotu, knotv, weight, W_xyi, dW_xyi_u, dW_xyi_v):\n",
    "    # % We pass in \n",
    "    # %   i = basis function number\n",
    "    # %   p = order or b-spline functions\n",
    "    # %   xi = the coordinate we are evaluating at (in parameter space)\n",
    "    # %   knot = the vector of knot values\n",
    "    # %   weights = a vector of all the 'weights' required for NURBS\n",
    "    N_ip_u, dN_ip_u = BsplineBasis(i, p, xi, knotu)\n",
    "    N_jq_v, dN_jq_v = BsplineBasis(j, q, yi, knotv)\n",
    "    numBasisFnsu = len(knotu) - 1 - p\n",
    "    numBasisFnsv = len(knotv) - 1 - q\n",
    "    if W_xyi == 0 and dW_xyi_u == 0 and dW_xyi_v == 0:\n",
    "        for ihat in range(numBasisFnsu):\n",
    "            for jhat in range(numBasisFnsv):\n",
    "                N_u, dN_u = BsplineBasis(ihat, p, xi, knotu)\n",
    "                N_v, dN_v = BsplineBasis(jhat, q, yi, knotv)\n",
    "                W_xyi = W_xyi + N_u * N_v * weight[(jhat - 1) * numBasisFnsu + ihat, 0]\n",
    "                dW_xyi_u = dW_xyi_u + N_v * dN_u * weight[(jhat - 1) * numBasisFnsu + ihat, 0]\n",
    "                dW_xyi_v = dW_xyi_v + N_u * dN_v * weight[(jhat - 1) * numBasisFnsu + ihat, 0]\n",
    "    R = N_ip_u * N_jq_v * weight[(j - 1) * numBasisFnsu + i, 0] / W_xyi\n",
    "    dR_u = weight[(j - 1) * numBasisFnsu + i, 0] * (W_xyi * dN_ip_u * N_jq_v - dW_xyi_u * N_ip_u * N_jq_v) / (W_xyi ^ 2)\n",
    "    dR_v = weight[(j - 1) * numBasisFnsu + i, 0] * (W_xyi * dN_jq_v * N_ip_u - dW_xyi_v * N_ip_u * N_jq_v) / (W_xyi ^ 2)\n",
    "\n",
    "    return R, dR_u, dR_v, W_xyi, dW_xyi_u, dW_xyi_v\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c0e68f9b",
   "metadata": {},
   "outputs": [],
   "source": [
    "def isod(No_nodeEl,xc,yc,p,q,knotVecu,knotVecv,w,numm,nd):\n",
    "    dRu = np.zeros(No_nodeEl,1)\n",
    "    dRv = np.zeros(No_nodeEl,1)\n",
    "    R  =  np.zeros(No_nodeEl,1)\n",
    "    W_xyi =0\n",
    "    dW_xyi_u=0\n",
    "    dW_xyi_v=0\n",
    "    for i in range(No_nodeEl):\n",
    "        jj,ii  = np.where(numm==nd[i])\n",
    "        R[i, 0], dRu[i, 0],dRv[i, 0], W_xyi, dW_xyi_u, dW_xyi_v = NURBSbasis_matlab( ii,jj, p,q ,xc,yc, knotVecu,knotVecv,w,W_xyi,dW_xyi_u,dW_xyi_v)\n",
    "    return dRu,dRv,R \n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "55b7fa45",
   "metadata": {},
   "outputs": [],
   "source": [
    "def isoJacobi1(No_nodeEl,dRu,dRv,xcoord,ycoord):\n",
    "    jacob1 = np.zeros(2, 2)\n",
    "    for i in range(No_nodeEl):\n",
    "        jacob1[0, 0] = jacob1[0, 0] + dRu[i] * xcoord[i]\n",
    "        jacob1[0, 1] = jacob1[0, 1] + dRu[i] * ycoord[i]\n",
    "        jacob1[1, 0] = jacob1[1, 0] + dRv[i] * xcoord[i]\n",
    "        jacob1[1, 1] = jacob1[1, 1] + dRv[i] * ycoord[i]\n",
    "    return jacob1\n",
    "\n",
    "def isoJacobi2(isonum):\n",
    "    jacob2 = np.zeros(2, 2)\n",
    "    jacob2[0, 0] = (isonum[1, 0] - isonum[0, 0]) / 2\n",
    "    jacob2[1, 1] = (isonum[3, 1] - isonum[1, 1]) / 2\n",
    "    return jacob2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7aa22f48",
   "metadata": {},
   "outputs": [],
   "source": [
    "def center(No_element, knotVecu, knotVecv):\n",
    "    element_center = np.zeros(No_element, 2)\n",
    "    for iel in range(No_element):\n",
    "        isonum = ISOnum(knotVecu, knotVecv, iel)\n",
    "        sum1 = 0\n",
    "        sum2 = 0\n",
    "        for point1 in range(4):\n",
    "            sum1 = sum1 + isonum[point1, 0]\n",
    "            sum2 = sum2 + isonum[point1, 1]\n",
    "\n",
    "        element_center[iel, 0] = sum1 / 4\n",
    "        element_center[iel, 1] = sum2 / 4\n",
    "\n",
    "    return element_center\n",
    "\n",
    "\n",
    "def isoelement_area(controlPts, knotVecu, p, knotVecv, q):\n",
    "    w = controlPts[:, 2]\n",
    "    nglx = 2\n",
    "    ngly = 2  # 2x2 Gauss quadrature\n",
    "    [element, elenode, numm] = ISOnumtotal(knotVecu, p, knotVecv, q, controlPts)  # 单元及控制点的编号\n",
    "    No_element = element.shape[0]  # 单元数\n",
    "    point2, weight2 = GaussPoint2(nglx, ngly)  # sampling points & weights\n",
    "    No_nodeEl = element.shape[1]  # 每个单元的控制点数\n",
    "    element_area = np.zeros(No_element, 2)\n",
    "    nd = np.zeros(No_nodeEl)\n",
    "    xcoord = np.zeros_like(nd)\n",
    "    ycoord = np.zeros_like(nd)\n",
    "    for iel in range(int(No_nodeEl)):\n",
    "        a = 0\n",
    "        for i in range(No_nodeEl):\n",
    "            nd[i] = element[iel, i]  # extract connected node for (iel)-th element\n",
    "            xcoord[i] = elenode[nd[i], 0]  # % extract x value of the node\n",
    "            ycoord[i] = elenode[nd[i], 1]  # % extract y value of the node\n",
    "        for intx in range(nglx):\n",
    "            x = point2[intx, 0]  # 高斯点横坐标\n",
    "            wtx = weight2[intx, 0]  # weight in x-axis\n",
    "            for inty in range(ngly):\n",
    "                y = point2[inty, 0]  # 高斯点纵坐标\n",
    "                wty = weight2[inty, 1]  # weight in y-axis\n",
    "                isonum = ISOnum(knotVecu, knotVecv, iel)  # 参数空间内的节点坐标\n",
    "                xc = (isonum[1, 0] - isonum[0, 0]) * x / 2 + (isonum[1, 0] + isonum[0, 0]) / 2  # 参数空间高斯点横坐标\n",
    "                yc = (isonum[3, 1] - isonum[1, 1]) * y / 2 + (isonum[3, 1] + isonum[1, 1]) / 2  # 参数空间高斯点纵坐标 \n",
    "                [dRu, dRv] = isod(No_nodeEl, xc, yc, p, q, knotVecu, knotVecv, w, numm, nd)  # 对参数空间求导\n",
    "                jacob1 = isoJacobi1(No_nodeEl, dRu, dRv, xcoord, ycoord)  # 物理空间到参数空间映射\n",
    "                jacob2 = isoJacobi2(isonum)  # 参数空间到高斯点母空间映射\n",
    "                detjacob1 = np.linalg.det(jacob1)  # determinant of Jacobian\n",
    "                invjacob1 = np.linalg.inv(jacob1)  # inverse of Jacobian matrix\n",
    "                detjacob2 = np.linalg.det(jacob2)  # determinant of Jacobian\n",
    "                a = a + wtx * wty * detjacob1 * detjacob2  # compute element matrix\n",
    "        element_area[iel, 0] = iel\n",
    "        element_area[iel, 1] = a\n",
    "    return element_area"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "67084ed5",
   "metadata": {},
   "outputs": [],
   "source": [
    "def ISO_shape(El_dof, No_nodeEl, element,\n",
    "              knotVecu, knotVecv, p, q, controlPts, numm, No_element, element_center):\n",
    "    nglx = 2\n",
    "    ngly = 2\n",
    "    w = controlPts[:, 2]\n",
    "    point2, _ = GaussPoint2(nglx, ngly)  # sampling points & weights\n",
    "    index = np.zeros(No_element, El_dof)\n",
    "    RI = np.zeros(No_element, nglx * ngly)\n",
    "    dRU = np.zeros(No_element, nglx * ngly)\n",
    "    dRV = np.zeros(No_element, nglx * ngly)\n",
    "    Re = np.zeros(No_element, 1)\n",
    "\n",
    "    nd = np.zeros(No_nodeEl)\n",
    "    for iel in range(No_element):  # loop for the total number of elements\n",
    "        for i in range(No_nodeEl):\n",
    "            nd[i] = element(iel, i)  # extract connected node for (iel)-th element\n",
    "        for intx in range(nglx):\n",
    "            x = point2(intx, 1)\n",
    "            for inty in range(ngly):\n",
    "                y = point2(inty, 2)\n",
    "                isonum = ISOnum(knotVecu, knotVecv, iel)\n",
    "                xc = (isonum(2, 1) - isonum(1, 1)) * x / 2 + (isonum(2, 1) + isonum(1, 1)) / 2\n",
    "                yc = (isonum(4, 2) - isonum(2, 2)) * y / 2 + (isonum(4, 2) + isonum(2, 2)) / 2\n",
    "                dRu, dRv, R1 = isod(No_nodeEl, xc, yc, p, q, knotVecu, knotVecv, w, numm, nd)\n",
    "                RI[iel, (intx - 1) * ngly + inty] = R1\n",
    "                dRU[iel, (intx - 1) * ngly + inty] = dRu\n",
    "                dRV[iel, (intx - 1) * ngly + inty] = dRv\n",
    "        _, _, R2 = isod(No_nodeEl, element_center(iel, 1), element_center(iel, 2), p, q, knotVecu, knotVecv, w, numm, nd)\n",
    "        Re[iel, 1] = R2\n",
    "    return RI, Re, dRU, dRV"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ad7d8699",
   "metadata": {},
   "outputs": [],
   "source": [
    "def ISO_Materialiso(v12,E1,Bt,G12,ceta):\n",
    "    Q = np.zeros(3, 3)\n",
    "    E2 = E1 / Bt\n",
    "    v21 = v12 * E2 / E1\n",
    "    Q[0, 0] = E1 / (1 - v12 * v21)\n",
    "    Q[0, 1] = v12 * E2 / (1 - v12 * v21)\n",
    "    Q[1, 0] = v21 * E1 / (1 - v12 * v21)\n",
    "    Q[1, 1] = E2 / (1 - v12 * v21)\n",
    "    Q[2, 2] = G12\n",
    "    cosd = lambda xceta: np.cos(xceta * np.pi / 180)\n",
    "    sind = lambda xceta: np.sind(xceta * np.pi / 180)\n",
    "    transform = np.array([\n",
    "        [cosd(ceta)*cosd(ceta), sind(ceta)*sind(ceta),-2*sind(ceta)*cosd(ceta)],\n",
    "        [sind(ceta)*sind(ceta), cosd(ceta)*cosd(ceta), 2*sind(ceta)*cosd(ceta)],\n",
    "        [sind(ceta)*cosd(ceta), -sind(ceta)*cosd(ceta), cosd(ceta)*cosd(ceta)-sind(ceta)*sind(ceta)]\n",
    "    ])\n",
    "    Materialiso = transform * Q * transform.T\n",
    "    return Materialiso"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9c42d846",
   "metadata": {},
   "outputs": [],
   "source": [
    "def ISO_matmtrx(ceta, lamuda, HT):\n",
    "    lamuda1 = lamuda\n",
    "    lamuda2 = lamuda / HT\n",
    "    cosd = lambda xceta: np.cos(xceta * np.pi / 180)\n",
    "    sind = lambda xceta: np.sind(xceta * np.pi / 180)\n",
    "    tran = np.array([[cosd(ceta), sind(ceta)],\n",
    "                     [-sind(ceta), cosd(ceta)]])\n",
    "    la = np.array([[lamuda1, 0],\n",
    "                   [0, lamuda2]])\n",
    "    matmtrx = tran.T * la * tran\n",
    "    return matmtrx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "310da768",
   "metadata": {},
   "outputs": [],
   "source": [
    "def isoDeriv2(No_nodeEl, dRu, dRv, invjacob1):\n",
    "    # --------------------------------------------------------------------------\n",
    "    #  Purpose:\n",
    "    #     determine derivatives of 2-D isoparametric shape functions with\n",
    "    #     respect to physical coordinate system\n",
    "    #  Variable Description:\n",
    "    #     dhdx - derivative of shape function w.r.t. physical coordinate x\n",
    "    #     dhdy - derivative of shape function w.r.t. physical coordinate y\n",
    "    #     No_nodeEl - number of nodes per element\n",
    "    #     dhdr - derivative of shape functions w.r.t. natural coordinate r\n",
    "    #     dhds - derivative of shape functions w.r.t. natural coordinate s\n",
    "    #     invjacob - inverse of 2-D Jacobian matrix\n",
    "    # --------------------------------------------------------------------------\n",
    "    dRx = np.zeros(No_nodeEl)\n",
    "    dRy = np.zeros(No_nodeEl)\n",
    "    for i in range(No_nodeEl):\n",
    "        dRx[i] = invjacob1[0, 0] * dRu[i] + invjacob1[0, 1] * dRv[i]\n",
    "        dRy[i] = invjacob1[1, 0] * dRu[i] + invjacob1(1, 1) * dRv[i]\n",
    "    return dRx, dRy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "55ece7b9",
   "metadata": {},
   "outputs": [],
   "source": [
    "def isoKine2D_F(No_nodeEl,dRx,dRy):\n",
    "    # --------------------------------------------------------------------------\n",
    "    #  Purpose:\n",
    "    #     determine the kinematic equation between strains and displacements\n",
    "    #     for two-dimensional solids\n",
    "    #  Variable Description:\n",
    "    #     No_nodeEl - number of nodes per element\n",
    "    #     dhdx - derivatives of shape functions with respect to x\n",
    "    #     dhdy - derivatives of shape functions with respect to y\n",
    "    # --------------------------------------------------------------------------\n",
    "    kinmtx2 = np.zeros(3, 2*No_nodeEl)\n",
    "    for i in range(No_nodeEl):\n",
    "        i1 = i * 2\n",
    "        i2 = i1+1\n",
    "        kinmtx2[0, i1] = dRx[i]\n",
    "        kinmtx2[1, i2] = dRy[i]\n",
    "        kinmtx2[2, i1] = dRy[i]\n",
    "        kinmtx2[2, i2] = dRx[i]\n",
    "    return kinmtx2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "efda3c13",
   "metadata": {},
   "outputs": [],
   "source": [
    "def isoEldof(nd, No_nodeEl,No_dof):\n",
    "    # --------------------------------------------------------------------------\n",
    "    #  Purpose:\n",
    "    #     Compute system dofs associated with each element\n",
    "    #  Variable Description:\n",
    "    #     index - system dof vector associated with element \"iel\"\n",
    "    #     nd - node connectivity for the (iel)-th element\n",
    "    #     No_nodeEl - number of nodes per element\n",
    "    #     No_dof - number of dofs per node\n",
    "    # --------------------------------------------------------------------------\n",
    "    index = np.zeros(No_nodeEl * No_dof)\n",
    "    k = 0\n",
    "    for i in range(No_nodeEl):\n",
    "        start = (nd[i] - 1) * No_dof\n",
    "        for j in range(No_dof):\n",
    "            k = k + 1\n",
    "            index[k] = start + j\n",
    "    return index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "26604088",
   "metadata": {},
   "outputs": [],
   "source": [
    "def ISO_lk_F(El_dof,No_nodeEl,element,elenode, No_element,dRU,dRV,knotVecu,knotVecv,No_dof,Bt,ceta):\n",
    "    E1 = 0.75e11 # elastic modulus\n",
    "    v12 = 0.3  #  Poisson's ratio\n",
    "    G12 = E1 / (2 * (1 + v12))\n",
    "    nglx = 2\n",
    "    ngly = 2\n",
    "    ko = np.zeros(No_element, nglx * ngly)\n",
    "    _, weight2 = GaussPoint2(nglx, ngly) #  sampling points & weights\n",
    "    index = np.zeros(No_element, El_dof)\n",
    "    matmtrx = ISO_Materialiso(v12, E1, Bt, G12, ceta)\n",
    "    for iel in range(No_element):  # loop for the total number of elements\n",
    "        k=np.zeros(El_dof,El_dof)   # initialization of element matrix to zero\n",
    "        nd = np.zeros(No_nodeEl)\n",
    "        xcoord = np.zeros_like(nd)\n",
    "        ycoord = np.zeros_like(nd)\n",
    "        for i in range(No_nodeEl):\n",
    "            nd[i] = element[iel, i]  # extract connected node for (iel)-th element\n",
    "            xcoord[i] = elenode[nd[i], 1]    # extract x value of the node\n",
    "            ycoord[i] = elenode[nd[i], 2]    # extract y value of the node\n",
    "        for intx in range(nglx):\n",
    "            wtx = weight2[intx, 0]      # weight in x-axis\n",
    "            for inty in range(ngly):\n",
    "                wty = weight2[inty, 1]   # weight in y-axis\n",
    "                dRu = dRU[iel, (intx - 1) * ngly + inty]\n",
    "                dRv = dRV[iel, (intx - 1) * ngly + inty]\n",
    "                isonum = ISOnum(knotVecu, knotVecv, iel)\n",
    "                jacob1 = isoJacobi1(No_nodeEl, dRu, dRv, xcoord, ycoord)\n",
    "                jacob2 = isoJacobi2(isonum)\n",
    "                detjacob1 = np.linalg.det(jacob1)  # determinant of Jacobian\n",
    "                invjacob1 = np.linalg.inv(jacob1)  # inverse of Jacobian matrix\n",
    "                detjacob2 = np.linalg.det(jacob2)  # determinant of Jacobian\n",
    "                dRx, dRy = isoDeriv2(No_nodeEl, dRu, dRv, invjacob1)\n",
    "                kinmtx2 = isoKine2D_F(No_nodeEl, dRx, dRy)\n",
    "                k = kinmtx2.T * matmtrx * kinmtx2 * wtx * wty * detjacob1 * detjacob2\n",
    "                ko[iel, (intx - 1) * ngly + inty] = k\n",
    "        index[iel,:]=isoEldof(nd, No_nodeEl, No_dof)\n",
    "    return ko,index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f64d00d3",
   "metadata": {},
   "outputs": [],
   "source": [
    "def isofbound(No_nodeEl,t,iel,num,knotVecu,knotVecv,element,elenode,numm,p,q,w):\n",
    "    point, weight = GaussPoint1(2)\n",
    "    f = np.zeros(No_nodeEl, 1)\n",
    "    isonum = ISOnum(knotVecu, knotVecv, iel)\n",
    "    nd = np.zeros(No_nodeEl)\n",
    "    xcoord = np.zeros_like(nd)\n",
    "    ycoord = np.zeros_like(nd)\n",
    "    for i in range(No_nodeEl):\n",
    "        nd[i]=element[iel,i]  # extract connected node for (iel)-th element\n",
    "        xcoord[i] = elenode[nd[i], 1]  # extract x value of the node\n",
    "        ycoord[i] = elenode[nd[i], 2]   # extract y value of the node\n",
    "    for int in range(2):\n",
    "        a = point[int]\n",
    "        wt = weight[int]\n",
    "        if num == 1:\n",
    "            y = a\n",
    "            x = -1\n",
    "        elif num == 2:\n",
    "            y = a\n",
    "            x = 1\n",
    "        elif num == 3:\n",
    "            y = -1\n",
    "            x = a\n",
    "        elif num == 4:\n",
    "            y = 1\n",
    "            x = a\n",
    "        xc = (isonum(2, 1) - isonum(1, 1)) * x / 2 + (isonum(2, 1) + isonum(1, 1)) / 2\n",
    "        yc = (isonum(4, 2) - isonum(2, 2)) * y / 2 + (isonum(4, 2) + isonum(2, 2)) / 2\n",
    "        R = np.zeros(1, No_nodeEl)\n",
    "        dRu = np.zeros(No_nodeEl, 1)\n",
    "        dRv = np.zeros(No_nodeEl, 1)\n",
    "        jac1 = 0\n",
    "        W_xyi = 0\n",
    "        dW_xyi_u = 0\n",
    "        dW_xyi_v = 0\n",
    "        jacx1 = 0\n",
    "        jacy1 = 0\n",
    "        for i in range(No_nodeEl):\n",
    "            jj, ii = np.where(numm == nd[i])\n",
    "            R[0, i], dRu[i, 0], dRv[i, 0], W_xyi, dW_xyi_u, dW_xyi_v = NURBSbasis_matlab(ii, jj, p, q, xc, yc, knotVecu,\n",
    "                                                                                         knotVecv, w, W_xyi, dW_xyi_u,\n",
    "                                                                                         dW_xyi_v)\n",
    "            if num == 1 or num == 2:\n",
    "                jacx1 = dRv[i, 0] * xcoord[i] + jacx1\n",
    "                jacy1 = dRv[i, 0] * ycoord[i] + jacy1\n",
    "            elif num == 3 or num == 4:\n",
    "                jacx1 = dRu[i, 0] * xcoord[i] + jacx1\n",
    "                jacy1 = dRu[i, 0] * ycoord[i] + jacy1\n",
    "        jac1 = np.sqrt(jacx1 ** 2 + jacy1 ** 2)\n",
    "        if num == 1 or num == 2:\n",
    "            jac2 = abs(isonum[3, 1] - isonum[1, 1]) / 2\n",
    "        elif num == 3 or num == 4:\n",
    "            jac2 = abs(isonum[1, 0] - isonum[0, 0]) / 2\n",
    "        f = f + R.T*abs(jac1)*wt*jac2*t\n",
    "    return f\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4a0eab86",
   "metadata": {},
   "outputs": [],
   "source": [
    "def ISO_broundary_F(knotVecu, knotVecv, q, p, element,\n",
    "                    elenode, numm, Sys_dof, controlPts, No_nodeEl):\n",
    "    uniqueKnots_u = np.unique(knotVecu)\n",
    "    uniqueKnots_v = np.unique(knotVecv)\n",
    "    n = uniqueKnots_u.shape[1] - 1\n",
    "    m = uniqueKnots_v.shape[1] - 1\n",
    "    nn = len(knotVecu) - 1 - p\n",
    "    mm = len(knotVecv) - 1 - q\n",
    "    No_element = element.shape[0]\n",
    "    kk = np.zeros(Sys_dof, Sys_dof)\n",
    "    ff = np.zeros(Sys_dof, 1)\n",
    "    U = np.zeros(Sys_dof, 1)\n",
    "    w = controlPts[:, 2]\n",
    "    rou = np.zeros(No_element, 4)\n",
    "    A = np.zeros(1, 2)\n",
    "    PP = np.array([0.18, 0])\n",
    "    C = len(elenode)\n",
    "    for i in range(len(elenode)):\n",
    "        A[0, :] = elenode[i, 1: 3]\n",
    "        if A == PP:\n",
    "            C = i\n",
    "    k = 1\n",
    "    D = np.zeros(1, 2)\n",
    "    bcdof = []\n",
    "    for j in range(len(elenode)):\n",
    "        D[1, :] = elenode[j, 1: 3]\n",
    "        if D[0, 1] == 0.18:\n",
    "            # bcdof[k] = j * 2 - 1\n",
    "            # bcdof[k + 1] = j * 2\n",
    "            # k = k + 2\n",
    "            bcdof.extend([j * 2 - 1, j * 2])\n",
    "    bcdof = np.array(bcdof)\n",
    "    bcval = np.zeros(k - 1, 1)\n",
    "    bcdof = [2 * np.arange(nn), 2 * (nn * (mm - 1) + np.arange(nn * mm)) - 1]  # first four dofs are constrained\n",
    "    bcval = np.zeros(4 * nn, 1)   # whose described values are 0\n",
    "    bcdof = [2 * np.arange(nn), 2 * np.arange(nn) - 1, 2 * (nn * (mm - 1) + np.arange(nn * mm)) - 1,\n",
    "             2 * (nn * (mm - 1) + np.arange(nn * mm))]  # first four dofs are constrained\n",
    "    bcval = np.zeros(4 * nn, 1)\n",
    "    nd = np.zeros(No_nodeEl)\n",
    "    for iel in range(n, n * m / 2 + n, n):\n",
    "        f = isofbound(No_nodeEl, 1e3 / 20, iel, 2, knotVecu, knotVecv, element, elenode, numm, p, q, w)\n",
    "        for i in range(No_nodeEl):\n",
    "            nd[i] = element[iel, i]   # extract connected node for (iel)-th element\n",
    "            ff[2 * nd(i) - 1, 0] = f[i, 0] + ff[2 * nd[i] - 1, 0]\n",
    "    nd = np.zeros(No_nodeEl)\n",
    "    for iel in range(n * (m / 2+1), n * m+n, n):\n",
    "        f = isofbound(No_nodeEl, 1.5e3 / 20, iel, 2, knotVecu, knotVecv, element, elenode, numm, p, q, w)\n",
    "        for i in range(No_nodeEl):\n",
    "            nd[i] = element(iel, i)  # extract connected node for (iel)-th element\n",
    "            ff[2 * nd(i), 1] = f(i, 1) + ff[2 * nd(i), 1]\n",
    "    ff[2 * C, 0] = -60000\n",
    "    return ff, bcdof, bcval\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "96cf85bf",
   "metadata": {},
   "outputs": [],
   "source": [
    "def ISO_F(index,ko,RI,d,penal,element,El_dof,No_element,ff,bcdof,bcval,Sys_dof):\n",
    "    iK = np.reshape(np.kron(index, np.ones(El_dof, 1)).T, (El_dof*El_dof*No_element, 1))\n",
    "    jK = np.reshape(np.kron(index, np.ones(1, El_dof)).T, (El_dof*El_dof*No_element, 1))\n",
    "    kk = np.zeros(Sys_dof, Sys_dof)\n",
    "    U = np.zeros(Sys_dof, 1)\n",
    "    sK1 = np.zeros(El_dof * El_dof, No_element)\n",
    "    return U\n",
    "def ISO_DC_F(No_nodeEl, U, No_element, No_dof, No_nodeSys, element, d, penal, ko, RI, El_dof):\n",
    "    dc = np.zeros(No_nodeSys, 1)\n",
    "    c = 0\n",
    "\n",
    "    for iel in range(No_element):\n",
    "        nd = np.zeros(len(No_element))\n",
    "        for i in range(No_nodeEl):\n",
    "            nd[i] = element[iel, i]  # extract connected node for (iel)-th element\n",
    "        index = isoEldof(nd, No_nodeEl, No_dof)\n",
    "        Ue = U(index, 1)\n",
    "        KE = np.zeros(El_dof, El_dof)\n",
    "        for i in range(4):\n",
    "            k = ko[iel, i]\n",
    "            R = RI[iel, i]\n",
    "            roud = R.T * d[element[iel, :], 0]\n",
    "            roug = penal * roud ^ (penal - 1)\n",
    "            KE = KE + roud ** penal * k\n",
    "            for I in range(No_nodeEl):\n",
    "                dc[nd[I], 0] = dc[nd[I], 0] - roug * R(I, 0) * Ue.T * k * Ue  # compute element matrix\n",
    "        c = c + Ue.T * KE * Ue\n",
    "    return dc, c\n",
    "\n",
    "\n",
    "def ISO_DV(No_element, No_nodeEl, Re, element_area, No_nodeSys, element):\n",
    "    dv = np.zeros(No_nodeSys, 1)\n",
    "    for iel in range(No_element):\n",
    "        R = Re[iel, 0]\n",
    "        for I in range(No_nodeEl):\n",
    "            dv[element[iel, I], 0] = dv[element[iel, I], 1] + R(I, 1) * element_area(iel, 1)\n",
    "    return dv\n",
    "\n",
    "def ISOcheck(d, dc, scale, BBB, BBBCI):\n",
    "    if scale == 0:\n",
    "        dcn = dc\n",
    "    else:\n",
    "        Nnodes = 9\n",
    "        rmin = scale * BBB[Nnodes, :]\n",
    "        Rmin = np.kron(rmin, np.ones(Nnodes, 1))\n",
    "        fac = Rmin - BBB\n",
    "        fac[np.where(fac < 0)] = 0\n",
    "        facsum = np.sum(fac, 0)\n",
    "        a = d[BBBCI] * dc[BBBCI] * fac\n",
    "        dcn = sum(a, 0)\n",
    "        dcn = dcn / (d * facsum.T)\n",
    "    return dcn\n",
    "\n",
    "def ISO_OC(total_area,d,volfrac,dc,dv,element_area, Re,No_element,element):\n",
    "    l1 = 0\n",
    "    l2 = 100000\n",
    "    move = 0.2\n",
    "    dnew = 0\n",
    "    V = 0\n",
    "    while (l2 - l1 > 1e-4):\n",
    "        lmid = 0.5 * (l2 + l1)\n",
    "        dnew = max(0.001, max(d - move, min(1., min(d + move, d * np.sqrt(-(dc / lmid / dv))))))\n",
    "        V = 0\n",
    "        for iel in range(No_element):\n",
    "            R = Re[iel, 0]\n",
    "            roue = R.T * dnew[element[iel,:],1]\n",
    "            V = V + roue * element_area[iel, 1]\n",
    "        if V - volfrac * total_area > 0:\n",
    "            l1 = lmid\n",
    "        else:\n",
    "            l2 = lmid\n",
    "    return dnew, V"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d2d813c9",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e690a52b",
   "metadata": {},
   "outputs": [],
   "source": [
    "def TOP_F(index_f, ko_f, RI, d, penal, element, El_dof, No_nodeEl, No_element,\n",
    "          ff, Re, element_area, No_nodeSys, elenode, scale, total_area, volfrac,\n",
    "          No_dof, bcdof_f, bcval_f, Sys_dof, LOOP, BBB, BBBCI):\n",
    "    C = []\n",
    "    L = []\n",
    "    D = []\n",
    "    loop = 0\n",
    "    while loop < LOOP:\n",
    "        loop = loop + 1\n",
    "        dold = d\n",
    "        f = 0\n",
    "        # ISO analysis\n",
    "        U = ISO_F(index_f, ko_f, RI, d, penal, element, El_dof, No_element, ff, bcdof_f, bcval_f, Sys_dof)\n",
    "        dc, c = ISO_DC_F(No_nodeEl, U, No_element, No_dof, No_nodeSys, element, d, penal, ko_f, RI, El_dof)\n",
    "        dv = ISO_DV(No_element, No_nodeEl, Re, element_area, No_nodeSys, element)\n",
    "        dc = ISOcheck(d, dc, scale, BBB, BBBCI)\n",
    "        d, V = ISO_OC(total_area, d, volfrac, dc, dv, element_area, Re, No_element, element)\n",
    "        # print result\n",
    "        change = max(abs(d - dold))\n",
    "        print([f'It.: {loop} Obj.:  {c} Vol.: {V / (total_area)}  ch.: {change}'])\n",
    "        C[loop] = c\n",
    "        L[loop] = loop\n",
    "        D[loop, :] = d\n",
    "    x = elenode[:, 1]\n",
    "    y = elenode[:, 2]\n",
    "    z = d\n",
    "    # 画图\n",
    "    # .......\n",
    "    return C, D\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "177bdf37",
   "metadata": {},
   "outputs": [],
   "source": [
    "def isoKine2D_T(No_nodeEl, dRx, dRy):\n",
    "    # --------------------------------------------------------------------------\n",
    "    #   Purpose:\n",
    "    #      determine the kinematic equation between strains and displacements\n",
    "    #      for two-dimensional solids\n",
    "    #   Variable Description:\n",
    "    #      No_nodeEl - number of nodes per element\n",
    "    #      dhdx - derivatives of shape functions with respect to x\n",
    "    #      dhdy - derivatives of shape functions with respect to y\n",
    "    # --------------------------------------------------------------------------\n",
    "    kinmtx2 = np.zeros(2, No_nodeEl)\n",
    "    for i in range(No_nodeEl):\n",
    "        kinmtx2[0, i] = dRx[i]\n",
    "        kinmtx2[1, i] = dRy[i]\n",
    "    return kinmtx2\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8000532d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def isoAplyc2(kk, ff, bcdof, bcval):\n",
    "    # --------------------------------------------------------------------------\n",
    "    #   Purpose:\n",
    "    #      Apply constraints to matrix equation [kk]{x}={ff}\n",
    "    #   Variable Description:\n",
    "    #      kk - system matrix before applying constraints\n",
    "    #      ff - system vector before applying constraints\n",
    "    #      bcdof - a vector containging constrained d.o.f\n",
    "    #      bcval - a vector containing contained value\n",
    "    # --------------------------------------------------------------------------\n",
    "    #   n=length(bcdof);\n",
    "    #   sdof=size(kk);\n",
    "    #\n",
    "    #   for i=1:n\n",
    "    #      c=bcdof(i);\n",
    "    #      for j=1:sdof\n",
    "    #         kk(c,j)=0;\n",
    "    #      end\n",
    "    #      kk(c,c)=1;\n",
    "    #      ff(c)=bcval(i);\n",
    "    #   end\n",
    "    n = len(bcdof)\n",
    "    kk[bcdof, :] = 0\n",
    "    for i in range(n):\n",
    "        kk[bcdof[i], bcdof[i]] = 1\n",
    "    ff[bcdof] = bcval\n",
    "    return kk, ff\n",
    "\n",
    "def ISO_T(index,ko,k_b,RI,d,penal,element,No_nodeEl,No_element,ff,No_nodeSys,elenode):\n",
    "    iK = np.reshape(np.kron(index, np.ones(No_nodeEl, 1)).T, (No_nodeEl * No_nodeEl * No_element, 1))\n",
    "    jK = np.reshape(np.kron(index, np.ones(1, No_nodeEl)).T, (No_nodeEl * No_nodeEl * No_element, 1))\n",
    "    temperature = np.zeros(No_nodeSys, 1)  # system force vector\n",
    "    kk = np.zeros(No_nodeSys, No_nodeSys)  # system matrix\n",
    "    sK1 = np.zeros(No_nodeEl * No_nodeEl, No_element)\n",
    "    for iel in range(No_element):  # loop for the total number of elements\n",
    "        ke = np.zeros(No_nodeEl, No_nodeEl)\n",
    "        for i in range(4):\n",
    "            k = ko[iel, i]\n",
    "            R = RI[iel, i]\n",
    "            roud = R.T*d[element[iel,:], 0]\n",
    "            ke = ke + roud ** penal * k\n",
    "            # ke=ke+k;\n",
    "        ke = ke + k_b[iel, 1]\n",
    "        sK1[:, iel]=ke[:]\n",
    "    sK = np.reshape(sK1, (No_nodeEl * No_nodeEl * No_element, 1))\n",
    "    kk = np.zeros([iK, jK, sK])  # todo: sparse matrix\n",
    "    kk = (kk + kk.T)/2\n",
    "    PC = [0, 0.45]\n",
    "    A = np.zeros(1, 2)\n",
    "    C = np.zeros(1, 1)\n",
    "    bcdof = 923\n",
    "    bcval = np.zeros(1, 1)\n",
    "    kk, ff = isoAplyc2(kk, ff, bcdof, bcval)\n",
    "    # temperature = kk\\ff\n",
    "    temperature = np.linalg.solve(kk, ff)\n",
    "    return temperature\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a20f06d6",
   "metadata": {},
   "outputs": [],
   "source": [
    "def ISO_lk_T(No_nodeEl, element, elenode, knotVecu, knotVecv, p, q, controlPts, numm, No_nodeSys, dRU, dRV, RI, Ht,\n",
    "             ceta):\n",
    "    uniqueKnots_u = np.unique(knotVecu)\n",
    "    uniqueKnots_v = np.unique(knotVecv)\n",
    "    n = uniqueKnots_u.shape[1] - 1\n",
    "    m = uniqueKnots_v.shape[1] - 1\n",
    "    nn = len(knotVecu) - 1 - p\n",
    "    mm = len(knotVecv) - 1 - q\n",
    "    No_element = element.shape[0]\n",
    "    ff = np.zeros(No_nodeSys, 1)  # system force vector\n",
    "    w = controlPts[:, 2]\n",
    "    lamuda = 2\n",
    "    matmtrx = ISO_matmtrx(ceta, lamuda, Ht)\n",
    "    nglx = 2\n",
    "    ngly = 2\n",
    "    ko = np.zeros(No_element, nglx * ngly)\n",
    "    k_b = np.zeros(No_element, 1)\n",
    "    w = controlPts[:, 2]\n",
    "    point2, weight2 = GaussPoint2(nglx, ngly)  # sampling points & weights\n",
    "    index = np.zeros(No_element, No_nodeEl)\n",
    "    for iel in range(No_element):  # loop for the total number of elements\n",
    "        #  initialization of element matrix to zero\n",
    "        k_b[iel, 1] = 0\n",
    "        nd = np.zeros(No_nodeEl)\n",
    "        xcoord = np.zeros_like(nd)\n",
    "        ycoord = np.zeros_like(nd)\n",
    "        for i in range(No_nodeEl):\n",
    "            nd[i] = element[iel, i]  # extract connected node for (iel)-th element\n",
    "            xcoord[i] = elenode[nd[i], 2]  # extract x value of the node\n",
    "            ycoord[i] = elenode[nd[i], 3]  # extract y value of the node\n",
    "        k = np.zeros(No_nodeEl, No_nodeEl)\n",
    "        f = np.zeros(No_nodeEl, 1)\n",
    "\n",
    "        for intx in range(nglx):\n",
    "            x = point2[intx, 0]\n",
    "            wtx = weight2[intx, 1]  # weight in x-axis\n",
    "            for inty in range(ngly):\n",
    "                y = point2[inty, 2]\n",
    "                wty = weight2[inty, 2]  # weight in y-axis\n",
    "                isonum = ISOnum(knotVecu, knotVecv, iel)\n",
    "                dRu = dRU[iel, (intx - 1) * ngly + inty]\n",
    "                dRv = dRV[iel, (intx - 1) * ngly + inty]\n",
    "                R1 = RI[iel, (intx - 1) * ngly + inty]\n",
    "                jacob1 = isoJacobi1(No_nodeEl, dRu, dRv, xcoord, ycoord)\n",
    "                jacob2 = isoJacobi2(isonum)\n",
    "                detjacob1 = np.linalg.det(jacob1)  # determinant of Jacobian\n",
    "                invjacob1 = np.linalg.inv(jacob1)  # inverse of Jacobian matrix\n",
    "                detjacob2 = np.linalg.det(jacob2)  # determinant of Jacobian\n",
    "                [dRx, dRy] = isoDeriv2(No_nodeEl, dRu, dRv, invjacob1)\n",
    "                kinmtx2 = isoKine2D_T(No_nodeEl, dRx, dRy)\n",
    "                k = kinmtx2.T * matmtrx * kinmtx2 * wtx * wty * detjacob1 * detjacob2  # compute element matrix\n",
    "                ko[iel, (intx - 1) * ngly + inty] = k\n",
    "                f = f + R1 * 0.1 * wtx * wty * detjacob1 * detjacob2\n",
    "\n",
    "        ff[element[iel, :], 0] = ff[element[iel, :], 0] + f[:, 0]\n",
    "        index[iel, :] = nd\n",
    "\n",
    "    return ko, k_b, index, ff\n",
    "\n",
    "def ISO_DC_T(No_nodeEl, temperature, No_element, No_nodeSys, element, d, penal, ko, k_b, RI):\n",
    "    dc = np.zeros(No_nodeSys, 1)\n",
    "    c = 0\n",
    "    for iel in (No_element):\n",
    "        nd = np.zeros(No_nodeEl)\n",
    "        for i in range(No_nodeEl):\n",
    "            nd[i] = element[iel, i]  # extract connected node for (iel)-th element\n",
    "        index = nd[:]\n",
    "        Te = temperature[index, 1]\n",
    "        KE = np.zeros(No_nodeEl, No_nodeEl)\n",
    "        for i in range(4):\n",
    "            k = ko[iel, i]\n",
    "            R = RI[iel, i]\n",
    "            roud = R.T * d[element[iel, :], 0]\n",
    "            roug = -penal * roud ** (penal - 1)\n",
    "            KE = KE + roud ** penal * k\n",
    "            for I in range(No_nodeEl):\n",
    "                dc[nd[I], 0] = dc[nd[I], 0] + roug * R[I, 0] * Te.T * k * Te  # compute element matrix\n",
    "        KE = KE + k_b[iel, 1]\n",
    "        c = c + Te.T * KE * Te\n",
    "    return dc, c\n",
    "\n",
    "\n",
    "def TOP_T(index_t,ko_t,k_b_t,RI,d,penal,element,No_nodeEl,No_element, pp,Re,element_area,No_nodeSys,elenode,scale,total_area,volfrac,LOOP,BBB,BBBCI):\n",
    "    S = []\n",
    "    L = []\n",
    "    D = []\n",
    "    loop = 0\n",
    "    while loop < LOOP:\n",
    "        loop = loop + 1\n",
    "        dold = d\n",
    "        f = 0\n",
    "        #  ISO analysis\n",
    "        temperature = ISO_T(index_t, ko_t, k_b_t, RI, d, penal, element, No_nodeEl, No_element, pp, No_nodeSys, elenode)\n",
    "        #  objective function and sensitivity analysis\n",
    "        ds, s = ISO_DC_T(No_nodeEl, temperature, No_element, No_nodeSys, element, d, penal, ko_t, k_b_t, RI)\n",
    "        dv = ISO_DV(No_element, No_nodeEl, Re, element_area, No_nodeSys, element)\n",
    "        ds = ISOcheck(d, ds, scale, BBB, BBBCI)\n",
    "        d, V = ISO_OC(total_area, d, volfrac, ds, dv, element_area, Re, No_element, element)\n",
    "        #  print result\n",
    "        change = max(abs(d - dold))\n",
    "        print(f'It.: {loop}    Obj.: {s}  Vol.:{V / (total_area)}   ch.: {change}')\n",
    "        S[loop] = s\n",
    "        L[loop] = loop\n",
    "        D[loop,:]=d\n",
    "    x = elenode[:, 1]\n",
    "    y = elenode[:, 2]\n",
    "    z = d\n",
    "    # todo: 畫圖\n",
    "    \n",
    "\n",
    "    return S,D\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "41295eaf",
   "metadata": {},
   "outputs": [],
   "source": [
    "#  控制点及单元显示\n",
    "# element_show(controlPts,knotVecu,knotVecv,p,q);\n",
    "volfrac=0.4\n",
    "penal=3.0\n",
    "No_dof=2\n",
    "No_nodeSys=controlPts.shape[0]                # total number of nodes in system\n",
    "# [element,elenode,numm]=ISOnumtotal(knotVecu,p,knotVecv,q,controlPts);%单元及控制点的编号\n",
    "bsFnConn_u, _ = iso_numelement(knotVecu,p)\n",
    "bsFnConn_v, _ = iso_numelement(knotVecv,q)\n",
    "\n",
    "numBasisFns_u=len(knotVecu)-1-p\n",
    "numBasisFns_v=len(knotVecv)-1-q\n",
    "num=np.zeros(numBasisFns_v,numBasisFns_u)\n",
    "\n",
    "for i in range(numBasisFns_v):\n",
    "    for j in range(numBasisFns_u):\n",
    "        num[i,j]=numBasisFns_u*(i-1)+j\n",
    "\n",
    "numm=num\n",
    "\n",
    "element = np.zeros(bsFnConn_u.shape[0]*bsFnConn_v.shape[0],bsFnConn_u.shape[1]*bsFnConn_v.shape[1])\n",
    "\n",
    "for i in range(bsFnConn_u.shape[0]):\n",
    "    for j in range(bsFnConn_v.shape[0]):\n",
    "        element[(j-1)*bsFnConn_u.shape[0]+i,:] = np.sort(num(bsFnConn_v[j,:],bsFnConn_u[i,:]).reshape(1,\n",
    "            bsFnConn_u.shape[1]*bsFnConn_v.shape[1]))\n",
    "        \n",
    "\n",
    "\n",
    "elenode=np.array([np.arange(numBasisFns_v*numBasisFns_u), controlPts])\n",
    "#  控制点距离\n",
    "Ax = elenode[:,0] - elenode[:,0].T\n",
    "Ay = elenode[:,1] - elenode[:,1].T\n",
    "A = np.sqrt(Ax ** 2 + Ay ** 2)\n",
    "B = np.sort(A)\n",
    "I = np.argsort(A)\n",
    "\n",
    "BBB=B[:9,:]\n",
    "BBBCI=I[1:9,:]\n",
    "\n",
    "#   参量初始化 \n",
    "\n",
    "No_element= element.shape[0]  # 单元数\n",
    "No_nodeEl= element.shape[1]   # 每个单元的控制点数\n",
    "El_dof=No_nodeEl*No_dof       # degrees of freedom per element\n",
    "Sys_dof=No_nodeSys*No_dof    # total system dofs \n",
    "\n",
    "d[:No_nodeSys, 0]=volfrac  # 所有控制点的密度值初始化\n",
    "nd = np.zeros(No_nodeEl,1)  # 向量初始化\n",
    "\n",
    "element_center=center(No_element,knotVecu,knotVecv)  # 参数空间单元中心点坐标\n",
    "element_area =isoelement_area(controlPts,knotVecu,p,knotVecv,q)\n",
    "element_area=element_area[:,1]  # 单元面积\n",
    "total_area=sum(element_area)  # 单元总面积\n",
    "loop = 0\n",
    "change = 1\n",
    "RI,Re,dRU,dRV  = ISO_shape(El_dof,No_nodeEl,element, knotVecu,knotVecv,p,q,controlPts,numm,No_element,element_center)\n",
    "#   热\n",
    "Ht=1\n",
    "ceta=0\n",
    "LOOP=151\n",
    "scale=1.2\n",
    "ko_t,k_b_t,index_t,pp = ISO_lk_T(No_nodeEl,element,elenode, knotVecu,knotVecv,p,q,controlPts,numm,No_nodeSys,dRU,dRV,RI,Ht,ceta)\n",
    "S_heat,D_heat = TOP_T(index_t,ko_t,k_b_t,RI,d,penal,element,No_nodeEl,No_element, pp,Re,element_area,No_nodeSys,elenode,scale,total_area,volfrac,LOOP,BBB,BBBCI)\n",
    "\n",
    "#    力\n",
    "Bt=1\n",
    "ceta=0\n",
    "LOOP=125\n",
    "scale=1.4\n",
    "ko_f,index_f  =  ISO_lk_F(El_dof,No_nodeEl,element,elenode, No_element,dRU,dRV,knotVecu,knotVecv,No_dof,Bt,ceta)  # 各个单元各个高斯点的刚度矩阵，形函数以及自由度\n",
    "ff,bcdof_f,bcval_f=ISO_broundary_F(knotVecu,knotVecv,q,p,element,elenode,numm,Sys_dof,controlPts,No_nodeEl)  # 边界条件\n",
    "C_static,D_static=TOP_F(index_f,ko_f,RI,d,penal,element,El_dof,No_nodeEl,No_element,ff,Re,element_area,No_nodeSys,elenode,scale,total_area,volfrac,No_dof,bcdof_f,bcval_f,Sys_dof,LOOP,BBB,BBBCI)\n",
    "U=ISO_F(index_f,ko_f,RI,d,penal,element,El_dof,No_element,ff,bcdof_f,bcval_f,Sys_dof)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.3 (tags/v3.10.3:a342a49, Mar 16 2022, 13:07:40) [MSC v.1929 64 bit (AMD64)]"
  },
  "vscode": {
   "interpreter": {
    "hash": "6f21bb2d5c565dbc33d815445cee0d4e3f2f7951520fd8c0e3b4200672f41bfe"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
